sink("R log.txt", split = TRUE)

library(ggplot2)
library(foreign)
library(plyr)
library(multiwayvcov)
library(lmtest)
library(lsmeans)
library(broom)
library(DataCombine)
library(zoo)
library(reshape2)
library(dplyr)
library(data.table)
library(tidyr)
library(boot)
library(ltm)
library(psych)
library(MCMCpack)
library(strucchange)

gc()




##################Appendix Figure 5#################

data <- read.dta("REPLICATION DATA.dta")

sdeez_ie <- NULL
sdeez_pe <- NULL
sdeez_cw <- NULL
sdeez_jg <- NULL

range_ie <- NULL
range_pe <- NULL
range_cw <- NULL
range_jg <- NULL

for (i in 1970:2014){
  
  loopdata <- data[data$year==i,]
  
  sdeez_ie[i] <- sd(loopdata$sub_lib_issueequal)
  sdeez_pe[i] <- sd(loopdata$sub_lib_policyequal)
  sdeez_cw[i] <- sd(loopdata$pollib_01)
  sdeez_jg[i] <- sd(loopdata$jake_pollib_01)

  range_ie[i] <- max(loopdata$sub_lib_issueequal) - min(loopdata$sub_lib_issueequal) 
  range_pe[i] <- max(loopdata$sub_lib_policyequal) - min(loopdata$sub_lib_policyequal)
  range_cw[i] <- max(loopdata$pollib_01) - min(loopdata$pollib_01)
  range_jg[i] <- max(loopdata$jake_pollib_01) - min(loopdata$jake_pollib_01)

}

sdeez_ie <- as.data.frame(sdeez_ie)
sdeez_pe <- as.data.frame(sdeez_pe)
sdeez_cw <- as.data.frame(sdeez_cw)
sdeez_jg <- as.data.frame(sdeez_jg)

names(sdeez_ie)[1] <- "stdev"
names(sdeez_pe)[1] <- "stdev"
names(sdeez_cw)[1] <- "stdev"
names(sdeez_jg)[1] <- "stdev"

sdeez_ie$scale <- "Issue Equal"
sdeez_pe$scale <- "Policy Equal"
sdeez_cw$scale <- "SPL"
sdeez_jg$scale <- "Expanded SPL"

sdeez_ie$year[1970:2014] <- 1970:2014
sdeez_pe$year[1970:2014] <- 1970:2014
sdeez_cw$year[1970:2014] <- 1970:2014
sdeez_jg$year[1970:2014] <- 1970:2014

sdeez <- rbind(sdeez_ie,sdeez_pe,sdeez_cw,sdeez_jg)
sdeez <- sdeez[is.na(sdeez$year)==F,]

sdeez$scale <- factor(sdeez$scale,levels(factor(sdeez$scale))[c(4,1,2,3)])


range_ie <- as.data.frame(range_ie)
range_pe <- as.data.frame(range_pe)
range_cw <- as.data.frame(range_cw)
range_jg <- as.data.frame(range_jg)

names(range_ie)[1] <- "stdev"
names(range_pe)[1] <- "stdev"
names(range_cw)[1] <- "stdev"
names(range_jg)[1] <- "stdev"

range_ie$scale <- "Issue Equal"
range_pe$scale <- "Policy Equal"
range_cw$scale <- "SPL"
range_jg$scale <- "Expanded SPL"

range_ie$year[1970:2014] <- 1970:2014
range_pe$year[1970:2014] <- 1970:2014
range_cw$year[1970:2014] <- 1970:2014
range_jg$year[1970:2014] <- 1970:2014

range <- rbind(range_ie,range_pe,range_cw,range_jg)
range <- range[is.na(range$year)==F,]

range$scale <- factor(range$scale,levels(factor(range$scale))[c(4,1,2,3)])

sdeez$measure <- "St.Dev."
range$measure <- "Range"

sdeezrange <- rbind(range, sdeez)


pdf(file="Appendix Figure 5.pdf", h=5,w=10)

ggplot(sdeezrange[is.na(sdeezrange$measure)==F,], aes(x=year, y=stdev, shape=scale)) +
  geom_point(size=1.5) +
  geom_line(size=.25) +
  scale_shape_manual(values=c(1,16,17,3), name="Scale") +
  #scale_linetype_manual(values=c(4,2,3,1), name="Scale") +
  facet_wrap(~measure) +
  xlab("Year") +
  ylab("") +
  scale_y_continuous(limits=c(0,1)) +
  theme_bw()

dev.off()




##########Appendix Figure 6###########

uni_files <- list.files(pattern="uni_")

mfx_list <- list()

for (i in seq_along(uni_files)){
  
  temp <- read.csv(uni_files[i], colClasses="character")
  temp <- temp[grepl("0.year2000to2014", temp$X)|grepl("1.year2000to2014", temp$X),]
  if(grepl("_dem", uni_files[i])){temp$party <- "Democratic"}
  if(grepl("_gop", uni_files[i])){temp$party <- "Republican"}
  if(grepl("_spl", uni_files[i])){temp$scale <- "SPL"}
  if(grepl("_expandedspl", uni_files[i])){temp$scale <- "Expanded SPL"}
  if(grepl("_policyequal", uni_files[i])){temp$scale <- "Substantive by Policy"}
  if(grepl("_issueequal", uni_files[i])){temp$scale <- "Substantive by Issue Area"}
  
  mfx_list[[i]] <- temp
  
}

plotdata <- do.call(rbind, mfx_list)

plotdata$mfx <- as.numeric(plotdata$.)
plotdata$se <- as.numeric(plotdata$X.1)

plotdata$cilow <- plotdata$mfx - (1.96*plotdata$se)
plotdata$cihigh <- plotdata$mfx + (1.96*plotdata$se)

plotdata$year[grepl("0.year", plotdata$X)] <- "1970-1999"
plotdata$year[grepl("1.year", plotdata$X)] <- "2000-2014"

plotdata$scale <- sapply(strwrap(plotdata$scale, 30, simplify=FALSE), paste, collapse="\n" )

plotdata$scale <- factor(plotdata$scale, levels(factor(plotdata$scale))[c(2,1,4,3)])

pd <- position_dodge(0.2)

pdf(file="Appendix Figure 6.pdf", w=12)

ggplot(plotdata, 
       aes(x=year, y=mfx, color=party)) +
  geom_point(position=pd) +
  geom_errorbar(aes(ymin=cilow, ymax=cihigh), width=0.1, position=pd) +
  facet_wrap(~scale, nrow=1) +
  geom_hline(yintercept=0, linetype=2) +
  scale_color_manual(name="Party\nControl", values=c("grey","black")) +
  ylab("Marginal Effect of Party Control") +
  xlab("Era") +
  theme_bw()

dev.off()





######################Figure 1###################

fullscales <- names(data[,grepl("fscale", colnames(data))==T & grepl("lag", colnames(data))==F])

data$abortionfscale[data$year<1973] <- NA
data$abortionfscale <- (data$abortionfscale - min(data$abortionfscale, na.rm=T)) / 
  (max(data$abortionfscale, na.rm=T) - min(data$abortionfscale, na.rm=T))

scales <- reshape(data, 
                  varying = c(fullscales), 
                  v.names = "score",
                  timevar = "scalearea", 
                  times = c(fullscales), 
                  direction = "long")


scales$scalearea <- as.character(scales$scalearea)
scales <- scales[!is.na(scales$scalearea) & !is.na(scales$state) & !is.na(scales$score) & !is.na(scales$year),]
scales$scalearea <- factor(scales$scalearea)


scales$scalename[scales$scalearea=="abortionfscale"] <- "Abortion"
scales$scalename[scales$scalearea=="campaignfinancefscale"] <- "Campaign Finance"
scales$scalename[scales$scalearea=="civilrightsfscale"] <- "Civil Rights"
scales$scalename[scales$scalearea=="criminaljusticefscale"] <- "Criminal Justice"
scales$scalename[scales$scalearea=="drugsfscale"] <- "Marijuana"
scales$scalename[scales$scalearea=="educationfscale"] <- "Education"
scales$scalename[scales$scalearea=="environmentfscale"] <- "Environment"
scales$scalename[scales$scalearea=="gunsfscale"] <- "Guns"
scales$scalename[scales$scalearea=="healthwelfarefscale"] <- "Health/Welfare"
scales$scalename[scales$scalearea=="housingtransmiscfscale"] <- "Housing/Transpo"
scales$scalename[scales$scalearea=="immigrationfscale"] <- "Immigration"
scales$scalename[scales$scalearea=="laborfscale"] <- "Labor"
scales$scalename[scales$scalearea=="publiclaborfscale"] <- "Public Sector Labor"
scales$scalename[scales$scalearea=="lgbtfscale"] <- "LGBT"
scales$scalename[scales$scalearea=="taxfscale"] <- "Tax"
scales$scalename[scales$scalearea=="votingfscale"] <- "Voting"



pdf(file="Figure 1.pdf", w=12,h=12)

ggplot(scales[scales$year>=1970 & is.na(scales$scalearea)==F & is.na(scales$score)==F &
                (scales$scalename!="Abortion"|scales$year>=1973),], 
       aes(x=year, y=score, group=state)) +
  geom_line(color="grey80") +
  geom_smooth(data=subset(scales, year>=1970 & is.na(partycontrol)==F &
                            is.na(scalearea)==F), 
              aes(x=year, y=score, group=factor(partycontrol), linetype=factor(partycontrol)),
              se=T,
              size=1, method="loess", color="black") +
  facet_wrap(~scalename) +
  ylab("Score (Left-Right)") +
  scale_linetype_manual(name="Party\nControl", values=c(3,1,2), labels=c("Republican","Divided","Democratic")) +
  coord_cartesian(ylim=c(-0.1,1)) +
  xlab("Year") +
  theme_bw()

dev.off()




######################Figure 2###################
marginsfiles <- list.files(pattern = "margins_dem_")

temp <- NULL
marginslist <- list()
cutoffrowlow <- NULL
cutoffrowhigh <- NULL

for (i in 1:length(marginsfiles)){
  
  temptemp <- read.csv(marginsfiles[i])
  
  tempdem <- temptemp[c(6,7),]
  tempdem$party <- "Democratic"
  temp <- tempdem
  
  temp$scalearea <- gsub("margins_dem_", "", marginsfiles[i])
  temp$scalearea <- gsub(".csv", "", temp$scalearea)
  
  temp$year <- c("1970-1999","2000-2014")
  
  names(temp)[names(temp)=="."] <- "effect"
  names(temp)[names(temp)=="X.1"] <- "se"
  
  marginslist[[i]] <- temp
  
}


marginsdem <- do.call(rbind.data.frame, marginslist)

marginsfiles <- list.files(pattern = "margins_gop_")

temp <- NULL
marginslist <- list()
cutoffrowlow <- NULL
cutoffrowhigh <- NULL

for (i in 1:length(marginsfiles)){
  
  temptemp <- read.csv(marginsfiles[i])
  
  tempgop <- temptemp[c(6,7),]
  tempgop$party <- "Republican"
  temp <- tempgop
  
  temp$scalearea <- gsub("margins_gop_", "", marginsfiles[i])
  temp$scalearea <- gsub(".csv", "", temp$scalearea)
  
  temp$year <- c("1970-1999","2000-2014")
  
  names(temp)[names(temp)=="."] <- "effect"
  names(temp)[names(temp)=="X.1"] <- "se"
  
  marginslist[[i]] <- temp
  
}

marginsgop <- do.call(rbind.data.frame, marginslist)

margins <- rbind(marginsdem, marginsgop)

margins$effect <- as.numeric(as.character(margins$effect))
margins$se <- as.numeric(as.character(margins$se))

margins$cilow <- margins$effect-(1.96*margins$se)
margins$cihigh <- margins$effect+(1.96*margins$se)

margins$model[grepl("lagdv", margins$scalearea)] <- "Lag DV"
margins$model[grepl("2way", margins$scalearea)] <- "2-Way FE"
margins$model[grepl("firstdiff", margins$scalearea)] <- "First Diff DV"

margins$scalearea <- gsub("lagdv_","", margins$scalearea)
margins$scalearea <- gsub("2wayfe_","", margins$scalearea)
margins$scalearea <- gsub("firstdiff_","", margins$scalearea)

margins$scalename[margins$scalearea=="abortionfscale"] <- "Abortion"
margins$scalename[margins$scalearea=="campaignfinancefscale"] <- "Campaign Finance"
margins$scalename[margins$scalearea=="civilrightsfscale"] <- "Civil Rights"
margins$scalename[margins$scalearea=="criminaljusticefscale"] <- "Criminal Justice"
margins$scalename[margins$scalearea=="drugsfscale"] <- "Marijuana"
margins$scalename[margins$scalearea=="educationfscale"] <- "Education"
margins$scalename[margins$scalearea=="environmentfscale"] <- "Environment"
margins$scalename[margins$scalearea=="gunsfscale"] <- "Guns"
margins$scalename[margins$scalearea=="healthwelfarefscale"] <- "Health/Welfare"
margins$scalename[margins$scalearea=="housingtransmiscfscale"] <- "Housing/Transpo/Misc"
margins$scalename[margins$scalearea=="immigrationfscale"] <- "Immigration"
margins$scalename[margins$scalearea=="laborfscale"] <- "Labor"
margins$scalename[margins$scalearea=="publiclaborfscale"] <- "Public Sector Labor"
margins$scalename[margins$scalearea=="lgbtfscale"] <- "LGBT"
margins$scalename[margins$scalearea=="taxfscale"] <- "Tax"
margins$scalename[margins$scalearea=="votingfscale"] <- "Voting"

pd <- position_dodge(0.2)


pdf(file="Figure 2.pdf", w=10,h=10)

ggplot(margins[margins$model=="Lag DV",], 
       aes(x=year, y=effect, color=party)) +
  geom_point(position=pd, size=2) +
  geom_errorbar(aes(ymin=cilow, ymax=cihigh), width=0, position=pd) +
  #geom_smooth(span=0.5, se=F) +
  facet_wrap(~scalename) +
  geom_hline(yintercept=0, linetype=2) +
  scale_color_manual(name="Party\nControl", values=c("grey","black")) +
  coord_cartesian(ylim=c(-0.04,0.04)) +
  #scale_y_continuous(breaks=c(-.04,-.02,0,.02,.04)) +
  ylab("Marginal Effect of Party Control") +
  xlab("Era") +
  theme_bw()

dev.off()




##################Figure 3###################

#Plot (a)
data$uninsured_pct <- 100*data$uninsured_pct

pdf(file="Figure 3a.pdf")

ggplot(data[is.na(data$uninsured_pct)==F & is.na(data$partycontrol)==F,], 
       aes(x=year, y=uninsured_pct, group=state)) +
  geom_line(color="grey80") +
  geom_smooth(aes(group=factor(partycontrol),linetype=factor(partycontrol)), color="black") +
  scale_x_continuous(breaks=c(1990,2000,2010), limits=c(1987,2014)) +
  ylab("Percent Uninsured") +
  scale_linetype_manual(name="Party\nControl", values=c(3,1,2), labels=c("Republican","Divided","Democratic")) +
  xlab("Year") +
  theme_bw()

dev.off()


#Plot (b)
uninsured <- read.csv("Party Effect on Uninsured Rate.csv")

uninsured$model <- sapply(strwrap(uninsured$model, 20, simplify=FALSE), paste, collapse="\n" )
uninsured$model = factor(uninsured$model,levels(factor(uninsured$model))[c(2,3,1)])

pd <- position_dodge(0.2)

pdf(file="Figure 3b.pdf")

ggplot(uninsured, 
       aes(x=year, y=mfx, group=interaction(model, party), shape=model, color=party)) +
  geom_point(size=2, position=pd) +
  geom_errorbar(aes(ymin=cilow, ymax=cihigh), width=0, position=pd) +
  geom_hline(yintercept=0, linetype=2) +
  scale_color_manual(name="Party\nControl", values=c("grey","black")) +
  #coord_cartesian(ylim=c(-0.05,0.05)) +
  #scale_y_continuous(breaks=c(-.05,0,.05)) +
  ylab("Effect of Party Control on Uninsured Rate") +
  scale_shape_manual(values=c(16,17,15), name="Model") +
  xlab("Year") +
  theme_bw()

dev.off()



########Incarceration Effect Plot#######

#Plot (a)
pdf(file="Figure 4a.pdf")

ggplot(data[is.na(data$incarcerate)==F & is.na(data$partycontrol)==F,], 
       aes(x=year, y=incarcerate, group=state)) +
  geom_line(color="grey60") +
  geom_smooth(aes(group=factor(partycontrol),linetype=factor(partycontrol)), color="black") +
  ylab("Incarceration Rate") +
  scale_linetype_manual(name="Party\nControl", values=c(3,1,2), labels=c("Republican","Divided","Democratic")) +
  xlab("Year") +
  theme_bw()

dev.off()


#Plot (b)
prison <- read.csv("Party Effect on Incarceration.csv")

prison$model <- sapply(strwrap(prison$model, 20, simplify=FALSE), paste, collapse="\n" )

prison$model <- factor(prison$model,levels(factor(prison$model))[c(2,3,1)])

pd <- position_dodge(0.2)

pdf(file="Figure 4b.pdf")

ggplot(prison[prison$population=="All",], 
       aes(x=year, y=mfx, group=interaction(model, party), shape=model, color=party)) +
  geom_point(size=2, position=pd) +
  geom_errorbar(aes(ymin=cilow, ymax=cihigh), width=0, position=pd) +
  geom_hline(yintercept=0, linetype=2) +
  scale_color_manual(name="Party\nControl", values=c("grey","black")) +
  #coord_cartesian(ylim=c(-0.05,0.05)) +
  #scale_y_continuous(breaks=c(-.05,0,.05)) +
  ylab("Effect of Party Control on Incarceration Rate") +
  scale_shape_manual(values=c(16,17,15), name="Model") +
  xlab("Year") +
  theme_bw()

dev.off()



sink()